Polydisperse Adsorption: Pattern Formation Kinetics, Fractal Properties, and 

Transition to Order 
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We investigate the process of random sequential adsorption of polydisperse particles whose size 
distribution exhibits a power-law dependence in the small size limit, P{R) ~ R a_1 . We reveal 
a relation between pattern formation kinetics and structural properties of arising patterns. We 
propose a mean- field theory which provides a fair description for sufficiently small a. When a — > oo, 
highly ordered structures locally identical to the Apollonian packing are formed. We introduce a 
quantitative criterion of the regularity of the pattern formation process. When a ^> 1, a sharp 
transition from irregular to regular pattern formation regime is found to occur near the jamming 
coverage of standard random sequential adsorption with monodisperse size distribution. 

PACS numbers: 81.10.Aj, 02.50.-r, 05.40.+j, 61.43.-j 



I. INTRODUCTION 

Random sequential adsorption (RSA) is an irreversible 
process in which particles are adsorbed sequentially and 
without overlaps and deposited particles cannot diffuse 
or desorb from the substrate. The RSA model has been 
initially applied to reactions along polymer chains [0. 
More recently, RSA processes have found a variety of 
other applications from adhesion of colloidal particles and 
proteins onto substrates Q to chemisorption || and epi- 
taxial growth for reviews see j^]. 

In this paper, we study a simple generalization of RSA 
that creates a rich dynamic behavior and results in com- 
plex spatial patterns. Namely, we consider polydisperse 
random sequential adsorption (PRSA) processes. Ad- 
sorption of mixtures has been addressed in a very few 
studies ll")!]. If a mixture contains a small number of 
different sizes, geometric and kinetic characteristics are 
primarily determined by the smallest size. In many appli- 
cations, however, the size distribution is continuous and 
spreads over several decades [|lQ[ . Therefore, before the 
smallest size will finally win, an interesting intermediate 
asymptotics arises. To address this intermediate regime, 
we consider PRSA with the power-law distribution in the 
small size limit, P(R) ~ i? a_1 . We show fll|] that this 
PRSA gives rise to fractal patterns of dimension Df de- 
pending on the exponent a. Additionally, we measure 
the degree of order of the patterns formed by PRSA, 
and identify the local structure of the patterns arising 
in the limit a — > oo with Apollonian packing [l^ ]. The 
significance of Apollonian packing in surface deposition 
problems was also recognized in Ref. ]l3[ . 

This paper is organized as follows. In Sec. II, we 
introduce PRSA and present numerical results. In 
Sec. Ill, we develop scaling, exact, and mean- field ap- 
proaches to PRSA. Exact results are available for the 
one-dimensional (ID) PRSA, and they are used to check 



mean- field and scaling approaches. A relation between 
kinetics and geometry of arising patterns is also dis- 
cussed. In Sec. IV, we introduce a quantitative criterion 
of regularity of the pattern formation process and ana- 
lyze the ordering in PRSA processes. The last Sec. V 
contains a summary. 



II. POLYDISPERSE RSA 

In applying RSA to real processes, one should take 
into account that adsorbed particles are typically poly- 
disperse. The relevant example is adsorption of colloidal 
particles, that have a broad radii distribution. It is usu- 
ally described by the Schulz distribution [Q, which has 
a power-law dependence on the radius R for small R and 
an exponential tail for large R: 



Psz(R) 



a 



R a ~ 



T(a) 



exp 



(R) 



R 



(1) 



Here (R) is the average radius and T(x) the gamma func- 
tion. Note that the exponent a should be positive to obey 
the normalization requirement, J dRP(R) = 1. 

Only the small-size behavior of P(R) affects the most 
interesting long time characteristics since in this regime 
only small particles can be adsorbed. Thus instead of (|l|) 
we shall use a power-law size distribution with an upper 
cutoff (taken as the unit of length): 



P(R) 



aR ' 1 R<1; 
R> 1. 



(2) 



The patterns formed by PRSA are drastically differ- 
ent from traditional RSA patterns, since the coverage is 
complete for PRSA. The pore space of the patterns is a 
nontrivial fractal set. This is physically evident, and in 
one dimension it proves possible to determine the fractal 
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dimension Df(a) analytically . In higher dimensions, we 
have to resort to numerical treatment. 

Monte Carlo simulations of the PRSA model have been 
performed by implementing the following algorithm. A 
center of the new disc is chosen at random with a uniform 
probability density. The radius of the disc is generated 
according to the size-distribution of Eq. (Q). If this disc 
does not overlap any other discs already in placed, it 
is deposited. Otherwise, the attempt is discarded. The 
maximal coverage studied in simulations was ~ 0.9 and 
about 100 runs were performed for each a. Some of the 
generated patterns are shown in Fig. 1. Clearly, the char- 
acter of patterns changes when the exponent a increases 
from to oo. For the small a, the patterns look like a 
random set of little discs distributed uniformly over the 
plane, with larger discs randomly scattered in the "sea" 
of smaller ones. For large a, one recognizes a structure 
initially formed by large discs and then reproduced by 
smaller discs in the holes between the large ones. These 
properties of patterns follow from the small size behav- 
ior of radii distribution function P(R). Indeed when 
a — > the smallest particles primarily participate in 
the adsorption, and hence apparently random patterns 
emerge. When a — > oo, the particles of the largest size 
are deposited until the system reaches the jamming limit 
of monodisperse RSA. Then the next particle to arrive 
will be the one that fits the biggest hole. This continues, 
so in this second stage the process is deterministic and 
apparently ordered patterns emerge. 

On the length scales much smaller than the upper cut- 
off, patterns are self-similar (Fig. 2). This suggests the 
fractal nature of the arising patterns. We have measured 
the fractal dimension Df of the pore-space of the pat- 
terns using the standard approach (see e.g. ]l5| , p!6| ) based 
on the analysis of the radii distribution function for the 
adsorbed particles. Namely, denote by n(R)dR the num- 
ber density of adsorbed discs with radii from the interval 
(R,R + dR). Let e is an arbitrary lower cutoff radius. 
Then N(e) = J°° dRn(R) gives the number density of 
discs with radii greater than e. The power law behavior 
of N(e) at the e — > limit, 



N(e) 



-D f 



(3) 



is a signature that the pore space is a set of fractal di- 
mension D f . Numerically, we indeed observed this power 
law behavior. We also found that Df monotonously de- 
creases when a increases (see Fig. 3). 



III. THEORETICAL APPROACHES TO 
PATTERN FORMATION 

In the present section we employ scaling, mean-field 
and exact approaches to the process of pattern forma- 
tion in PRSA. The first approach is based on the scaling 
hypothesis and gives relations between structural and ki- 
netic characteristics of the system. 



A. Scaling Framework 

Let is the fraction of uncovered area and \I/(i?, t) 
the probability that a disc of radius R can be placed onto 
a substrate IO]. Clearly, 



*(0,t) = *(*). 
evolves according to the exact rate equation 
d<S> 



dt 



dRP(R)<t>(R,t)V d R d , 



(4) 



(5) 



where Vd denotes the volume of the d-dimensional unit 
sphere, Vd = 7r d / 2 /r(l + d/2). Assuming a scaling be- 
havior of \I/(i?, i), we write 



n^t) = s v {t)F[-^j 



(6) 



Here S(t) ~ t~ v is a typical gap between neighboring 
adsorbed particles and F(x) is a scaling function. The 
scaling description applies when t — > oo and R — > with 
R/S(t) finite. The existence of scaling is an assumption, 
which is supported by numerical evidence in 2D and by 
analytical results in ID ||. 

Eqs. (|) and (0) imply $(t) ~ t~ z with z = 6v. Sub- 
stituting then (|) into (||) gives 



aVd / dx x' 
Jo 



.at+d-1 



F{x). 



(7) 



Scaling suggests that self-similar fractal structures 
should arise. Indeed, computing the number density of 
the absorbed particles we get 

POO 

n(R) = / dtP(R)V(R,t) ~ J R«- 1 +(^- 1 )/^ (8) 
Jo 

Hence, the power-law dependence of Eq. (||) is recovered 
when we identify the fractal dimension with 



Df — d — z(d + a). 



(9) 



Eq. ([)]) provides a relation between the fractal dimension 
of the arising patterns and the kinetic exponent z. Sim- 
ilar qualitative behaviors were observed in other pattern 
formation models p6p . 

We should stress that scaling provides just a frame- 
work; for instance, it gives scaling relations among the 
exponents but it does not allow to compute the expo- 
nents. So one should use other approaches to get a com- 
plete exact or approximate description of PRSA. In one 
dimension, an exact description is indeed possible. In 
higher dimensions, even in the extreme case of a = oo 
the fractal dimension of the pore space remains unknown, 
so analytical description of PRSA is hardly possible. 
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PRSA. a=0.1 , D f = 1.94 




PRSA. o=5 , Df = l.56 




PRSA. a=50 , I^ = 1 - 33 



Fig. 1 Typical PRSA patterns for a = 0.1, 1, 5, 




PRSA. a= 10 , Zy = 1-50 




Apollonian packing 



, 50 and Apollonian packing. 
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Fig. 2 Self-similarity of the PRSA patterns. 



B. ID PRSA: Exact Results 

A detailed analysis of ID PRSA is given in Ref. ||. 
Here we sketch basic results which are necessary to de- 
termine the fractal dimension. Let C(x,t) be the con- 
centration of holes of length x at time t. All holes have 
the same "shape" in ID that significantly simplifies the 
problem. By definition, C(x,i) is related to *&(x,t) via 



poo 

/ dy(y 

J X 



x)C(y,t). 



(10) 



The kinetic equation for C(x, t) reads 
dC(x,t) 



at 



-C(x,t) 



/ 

Jo 



dz(x - z)P(z) 



+2 



roo ry-x 
/ dyC(y,t) / dzP(z). 

Jx JO 



(11) 



The first term on the right-hand side of Eq. ( |l l| ) gives 
the loss of holes of length x due to deposition of intervals 
of length z (with z < x)\ the second term describes the 
gain of holes of length x from larger holes. Substituting 
P(z) = P(R) given by Eq. (§) into Eq. (O) yields 



d_ 

dt 



a + 1 



C{x,t) 



dyC(y,t)(y-x) a . (12) 



Multiplying both sides of Eq. ( |l2|) by x@ and integrating 
over x gives the kinetic equation 



dt 



r(a + l)rQ3 + l) 1 
T(a + /? + 2) a + 



M a+f3+1 (13) 



for the moments of the hole-size distribution 



rOO 

M (t) = / dxx p C(x,t). (14) 
Jo 

Eq. ([To]) implies C(x,t) = -j^ ty{x,t). Combining this 
with the scaling ansatz of Eq. (H) we obtain 



C(x,t) =S"- Z (t) F (x/S(t)), 



(15) 



where = d 2 F/d^ 2 . Eq. (|T^) allows us to express 

the moments Mp(t) via the time-independent moments, 

M {t) = S 6+ P- 1 {t)m p , m = d£^F"(0. (16) 

JO 

Choose now the exponent (3 so that the numerical factor 
on the right-hand side of Eq. (N3) vanishes, i.e., 



r(a + l)r(/3 + 1) 1 



T(a + /5 + 2) 



(17) 



For such (3 = /3(a), Eq. (|l3|) implies that Mp(t) does not 
depend on t. Eq. ( [l6| ) therefore gives = 1 — /3, and then 
other exponents and the fractal dimension are found: 



1 



1 + a 



A simple analysis shows that (i) Eq. ( |17| ) has only one 
positive solution (3 = /3(a), (ii) (3 < 1 for all a > 0, 
and (iii) /? decreases when a increases. One can de- 
termine (3 = Df explicitly in some specific cases, e.g., 
Df = (vTr — 3)/2 for the uniform size distribution 
(a = l). 
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C. Mean-Field Theory of PRSA 

In higher dimensions we employ an approximate mean- 
field treatment. We shall use a mean-field descrip- 
tion of PRSA close to the one developed in ]l6[ for 
nucleation-and-growth processes. As usual, we shall ig- 
nore many-particles spatial correlations and account only 
two-particle ones. To find t), consider a circle of ra- 
dius R centered at the origin. Then one can write the fol- 
lowing estimate to this function in the mean- field spirit: 



exp< - 



t nR 

dr / dr Vt^ 
o Jo 



d-1 



> * p( ')^T (19> 



x - p {-jr*r* n **" i jC*^ ) w}- 

The former exponential factor in ( |l9| ) estimates the prob- 
ability that our circle is not covered up to time t by discs 
with centers fall inside this circle. The latter exponential 
factor guarantees that the free space inside our circle is 
not covered by the other discs which were adsorbed out- 
side the circle. We also denote by the surface area of 
unit sphere in d dimensions, fl^ = dVd — 2 / K d ^ 2 /T(d/2). 
The exponential factors in Eq. (|l9|) were derived by not- 
ing that 

dr dr dp ft d r d_1 P(p)*(p, r) 

gives the probability that a disc of radius belonging to the 
interval (p, p + dp) was adsorbed within the time interval 
(r, r + dr) in the spherical shell centered at the origin 
and confined by radii r and r + dr. The factor r) 
guarantees that the disc of radius p may be placed into 
the system. Such event prevents the adsorption of disc of 
radius R at the origin. The probability that such event 
has not happened, 

1 - dr dr dp n d r d - 1 P(p)^{p, r), 

may be re-written as 

exp {-dr dr dp^r^ 1 P(p)V(p,r)} . 

The probability that none of these events has happened 
up to time t is obtained by multiplying all these factors 
with 0<r<£, 0<r<oo, and < p < oo. However 
one should take into account that above expression gives 
a correct estimate only for factors with r = 0. For sub- 
sequent factors with r > 0, one should use t)/$(t) 
instead of r) since the preceding factors guarantee 
that the discs are placed on the uncovered space. Using 
<3>(0) = 1 and treating separately r < R and r > R, one 
arrives at Eq. (|l9|). Again we assume that in the scal- 
ing regime we can use Eq. (||) for the function ^(R,t). 
Taking into account that ^(0,t) = we get 



F{Rt v ) 
exp • 



(20) 



dr I dp n d P(p) {(R + p) d - p d ) F(prn) ■ 



Remarkably, the ansatz 

F(x) = exp(— A\x ■ 



A d x d ) 



(21) 



reduces the nonlinear integral equation (g0|) to a system 
of d algebraic equations for coefficients Aj. In particular, 
in ID we get v = l/(a + 1) and a closed equation for A±: 



2a 



1 



VOL 



f 

Jo 



dx x^e'^. 



(22) 



Eq. (||) is solved to find A x = [2T(a + 2)] . Thus on 
the mean-field level, the scaling function for ID PRSA 
is pure exponential, i.e. it is clearly different from the 
analytical solution. To obtain a quantitative difference 
we substitute F — e~ AlX with the above value of A\ into 
Eq. (0) to find z = Substituting this into Eq. (^) 

we obtain the fractal dimension: Dy F = 1 — a for a < 1 
and Dy F = for a > 1. Thus the mean- field theory 
is clearly wrong for a > 1, though in the limit a — > 



it becomes exact [0. For instance, £)e xact — D 



MF 



f ~ f 

(4f + 4 7 2 - 2 - 67)a 2 + G(a 3 ) (here 7 ^ 0.57721566 is 
the Euler constant). 



For 2D PRSA, the ansatz of Eq. (|2lj) yields 

pOO 

A ± = 2ira(a + 2) / dx x a e~ AlX - A2X \ 
Jo 



roo 

A 2 = -a(a + 2) / dxx a ~ x t 
2 Jo 



-A\x— A 2 x 



Solving these equations numerically, and then inserting 
F(x) = e - A ^ x ~ A ^ x into Eq. (]?]), one finds z and Df. 

In the small a limit, we perform a perturbation anal- 
ysis to find Q 

A\ = 2-kol + a\o? — 2ira — 17.3557a 2 , 

A 2 = 7T + a 2 a + a 3 a 2 = tt - 1.8339a + 12.6546a 2 , 

where we have omitted terms of order 0(a 3 ). Using these 
expressions for A\ and A 2l we get the kinetic exponent 



z = a/2 + a 4 a 2 = a/2 - 1.5428a 2 



(23) 



and for the fractal dimension 

Df = 2-a J r a 5 a 2 = 2 - a + 2.5856a 2 + • • • . (24) 

The mean- field predictions for Df(a) and z(a) are 
shown on Fig. 3 and Fig. 4. We see that the mean-field 
approach provides a fair description for small a. For 
a ~ 1, however, the spatial correlations become more 
and more important and the mean- field theory fails. 
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Fig. 3 The fractal dimension Df vs. a in two dimensions. 
Df of the scaling theory is obtained from Eq. (^|) with 
the "experimental" value of the kinetic exponent z. 
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Fig. 4 The kinetic exponent z vs. a in one and two di- 
mensions. 



D. The Limit a — > oo 

When a = oo, the PRSA process develops through 
two stages. The initial stage is just RSA of monodis- 
perse particles of radii R = i? max = 1. It continues until 
the jamming coverage, = 0.542, is reached. Then 
the late stage begins, where the next disc to place would 
be the one that fits the biggest hole. The dynamics in 
this late stage is thus deterministic and extremal. This 
deterministic procedure has been applied by Apollonius 
of Perga 200BC to fill the space between three "kiss- 
ing" discs iQ. In the present case, the procedure fills 
uncovered space obtained during the initial RSA stage. 
As the process develops, new discs will be placed more 
and more often into the curvilinear triangles confined by 
three kissing discs. Such curvilinear triangles are filled 
independently, so locally our packing should be identical 
to Apollonian packing. The fractal dimension quantifies 



local characteristics of the pattern, so we conclude that 
D f (oc) = D A . 

The fractal dimension of Apollonian packing is Da — 
1.3057 in two dimensions. Note that the fractal di- 
mension of the Apollonian parking, arguably the old- 
est known fractal, has not been computed analytically 
in higher dimensions, so the exact value of Df (oo) = Da 
remains unknown. The only trivial exception is the one- 
dimensional case where the holes remaining after the ini- 
tial stage are filled up completely during the determin- 
istic stage. The number density n(x) = C(x,t = oo) of 
adsorbed intervals of length x is p9|| 



n(x) 



Jo 



dt t exp 



-xt 



Jo 



dr- 



1 



(25) 



n(x) exhibits a weak integrable singularity in the small 
size limit, n(x) ~ ln(l/x), implying that N(e) is regular 
and therefore Df (oo) = |p0| . This agrees with previ- 
ous exact results, Eqs.(17) and (jig), which provide the 
asymptotic behavior Df(a) ~ when a — > oo. 



ln(a+2) 



IV. ORDERING IN THE RSA 

For a quantitative description of the ordering processes 
we introduce an entropy, Sn, characterizing the degree of 
order of Af-particle patterns. Let C denotes a particular 
Af-particle pattern and Pn(C) denotes the probability of 
that pattern. We can define the Shannon entropy [pi|,^2[: 



S N = - ^2p N (C) log 2 p N (C) 



(26) 



As it follows from Eq. (|26|), Sn = for a regular pat- 
tern, since only one configuration (which occurs with the 
probability pat = 1) contributes to the entropy. A closely 
related quantity, dSN/dN ~ Sn+i — Sn, gives the en- 
tropy production rate per particle and characterizes the 
regularity of the pattern formation process. We also in- 
troduce the conditional entropy, £^ +1 (C), characterizing 
patterns built by adding a disc to a given Af-particle pat- 
tern C: 

S* N+1 (C) = -J2p(R^\C)log 2 p(R,f\C). (27) 

R,r 

Here p(R, r\C) is the conditional probability to add a 
disk of radius R at point r to the particular pattern C(N) 
of N disks. The total probability of the N + 1-particle 
configuration, obtained from the pattern C by placing an 
additional disc, reads: 



p N +i(R, f , C) = p(R, r\C)p N (C). 
Sn+i can be written as 



(28) 



S N+ i = - Y, PN+i(R,r,C)log 2 p N+1 (R,r,C). (29) 

R,f,C 
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Using Eqs. (|26|)-(|29|) and the normalization condition, 
rP(Ri r\C) = 1, we finally arrive at the entropy pro- 
duction rate: 



dS 



N 



dN 



Sn+i — Sn 



c 



p N (C)S* N+1 (C). (30) 



Thus the entropy production rate is obtained by averag- 
ing the conditional entropy over all possible configura- 
tions. 
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Fig. 5 The entropy production rate vs. the uncovered 
area for different values of a. A sharp decay at « 0.55 
is clearly seen for large a. 

In spite of apparent simplicity of the definition of the 
entropy production rate, its analytical evaluation is a 
very challenging problem even in one dimension. In- 
deed, it requires knowledge of the multiparticle proba- 
bility distribution function. In the simplest case of the 
monodisperse RSA we do know the first two moments 
of p(C), the density || and the pair correlation func- 
tion , but much more detail information is needed for 
determination of the entropy production rate. Thus we 
investigated this quantity numerically for 2D system by 
means of Monte Carlo simulations. We implemented the 
following algorithm. A disc which is added to the pat- 
tern C(N) of N disks was treated as a point in a the 
"configurational" space (R,r) = (R,x,y). We divided 
the continuous configuration space into sufficiently small 
discrete cells and enumerated these cells. For computa- 
tions we used cells of linear size 0.01 and 5x5 fragment 
of the surface (with R ma x = 1, see (^J)); thus, the total 
number of cells was M = 2.5 x 10 T . Then the conditional 
probabilities, p(R,f\C(N)) = p i: that (N + l)-th disk 
comes to the i th cell, corresponding to (i?, r) in the given 
configuration, C(N), of N discs was calculated numeri- 
cally from 



AiRT 



(31) 



where Ai = if the disc corresponding to the i th cell 
overlaps with some disc in the pattern C(N); otherwise, 
Ai = 1. Eq. ( |3l| ) follows from the rate of adsorption, 
(||), and the normalization condition for the probability. 
The discrete conditional probabilities ( |3l| ) are used to 
calculate the conditional entropy S N+1 (C) via Eq. (|2^). 
Finally, the entropy production rate was obtained by av- 
eraging S N+1 (C) over C according to Eq. (|30|). The av- 
eraging was performed over 10 2 Monte Carlo runs, and 
the accuracy of the method was controlled through run- 
to-run deviations. 

To compare the entropy production rate for different 
values of a, we plot dSw/dN versus (see Fig. 5). The 
striking behavior of the entropy production rate is clearly 
seen in the large a limit: At the beginning of the pro- 
cess of pattern formation (i.e., when « 1), the entropy 
production rate decreases slowly similar to the small a 
case. Around ~ 0.55, however, a sharp decay occurs. 
The threshold value, ~ 0.55, is close to the jamming 
density, ^ = 0.542, of the RSA of identical discs. The 
transition demarcates the initial random stage when discs 
of maximal radii are deposited and the final deterministic 
stage when discs of maximal possible radius are inserted 
into maximal holes. Given the deterministic nature of the 
final stage, the entropy production rate should be equal 
to zero when a = oo and $ > ^oo. This is in agreement 
with our numerical findings for large a. 



V. SUMMARY 

In summary, we investigated the adsorption kinetics 
and spatial properties of the arising patterns in RSA 
of polydisperse particles whose size distribution has a 
power-law form in the small size limit. We developed 
a scaling approach and verified that it indeed applies 
by comparing with exact results in one dimension and 
numerical results in two dimensions. We found that 
arising patterns are self-similar fractals that appear to 
be completely random when deposited particles are pre- 
dominantly small; in the opposite limit highly ordered 
structures, locally isomorphic to Apollonian parking, are 
formed. The fractal dimension Df of the pore space is 
determined by the power-law exponent a of the particles 
size distribution. When a increases from to oo, Df de- 
creases from 2 to Da — 1.305, of the Apollonian packing. 
We introduced the entropy production rate as a quanti- 
tative measure of the regularity of arising patterns. We 
observed that for sufficiently small a, the entropy pro- 
duction rate smoothly decays as the coverage increases. 
In the complimentary case of a ^> 1, the entropy produc- 
tion rate displays a similar behavior for sufficiently small 
coverage, followed by a sharp decay to a very low en- 
tropy production rate. Physically, it reflects a two-stage 
nature of the pattern formation process in the large a 
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limit: The ordinary RSA, that goes until the jammed 
state is reached, leads to a random structure which is 
starting point for the second stage. During this late stage 
the deposition process is deterministic and extremal - at 
a given step, the largest hole is filled. 

One of us (PLK) acknowledges NSF grant DMR-9632059 
and ARO grant DAAH04-96-1-0114 for financial support. 
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